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ABSTRACT 


This thesis focuses on producing realistic time-space-position information 
(TSPI) data for moving targets as inputs for a weapon flyout model to compute 
the weapon Circular Error Probable (CEP) and the probability of damage. The 
velocity profile of the target is modeled based on the kinematic constraints for the 
type of vehicle and the type of path on which it is traveling. The discrete-time 
position profile can then be derived for the target and fed into the weapon flyout 
model to compute the associated measures of weapon effectiveness. The 
methodology is applied to both land vehicles traveling on specific roads, 
differentiated based on the road profiles, and small boats moving on the open 
sea to generate various sets of measures of weapon effectiveness for different 
scenarios. A program was also developed to enable a quick classification of the 
target’s expected route of advancement into one of the six road types. This 
allows the user to perform quick analysis on the target’s planned route of travel 
and produce measures of weapon effectiveness, which can then be readily 
utilized for swift decision making on the weapon options available. 
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I. INTRODUCTION 


A. MOTIVATION 

In modern warfare, with the rapid development of military technologies, 
armed forces around the world are increasingly using guided weapons to protect 
their valuable assets and to destroy targets of interest for strategic gains. Such 
weapons can be launched from various platforms, including aerial, surface, or 
even submerged ones. Even with the proliferation of the guided weapons, the 
ability of the weapons to hit moving targets remains a key area of interest. 
According to [1], the U.S. Navy is in the midst of equipping the submarine and 
the Tomahawk missile with enhanced seeker technology in order to better enable 
the weapon to take out mobile naval threats of the hostile forces. Meanwhile, [2] 
reveals that even the F-35, despite being the state-of-the-art fifth generation 
stealth fighter, is still looking for options to enhance its capability to hit moving 
targets. 

Depending on the type of target environment, and the target’s 
vulnerabilities, extensive weapon assessment is required to evaluate the likely 
weapon performance against such target, and the weapon delivery profile for the 
best possible outcome. To perform such assessment with models of higher 
fidelity, it is important to achieve a realistic representation of the target motion 
based on the target’s type of environment. 

This thesis focuses on the position estimation for both a fast-moving land 
vehicle along an expected route and a fast-moving boat on an open sea. The aim 
of the thesis is to develop a methodology to model an estimator for the likely path 
of travel for the land and sea targets of interest, and for the land target, to classify 
the target environment based on the road profile and vehicle characteristics; 
using this information enables a set of weapon performance estimates to be 
computed to facilitate quick decision making at the operational level. The goal of 
this thesis is to develop a realistic model characterizing the most probable paths 
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of interest for specific target types based on the target environment to provide 
useful inputs that can be used to generate weapon performance metrics that 
enable the choice of weapon delivery for best possible outcome. 

B. BACKGROUND 

Weaponeering is defined in [3] as the process of “determining the quantity 
of a specific type of weapon required to achieve a defined level of target damage, 
considering target vulnerability, weapon effects, munitions delivery error, damage 
criteria, probability of kill, weapon reliability, etc.” By determining the statistical 
estimate of the weapon’s effectiveness against specific target, the user can then 
designate the appropriate type and quantity of weapons to be used to achieve 
the level of damage required. For weaponeering to achieve a close estimate of 
the expected weapon performance, it is important for the parameters assumed to 
match those used during weapon deployment in combat. Otherwise, an 
inaccurate model will lead to a poor weaponeering assessment, and result in 
inefficient or ineffective use of the weapons. 

The area of interest is in making realistic weaponeering assessment 
against moving targets. Broadly, this can be classified into two scenario types, 
constrained and unconstrained. The first scenario type involves targets moving 
along a defined path, for example, a vehicle going down a particular stretch of 
road or a train on a railroad track, while the latter scenario involves targets 
moving in an open environment, for example, a tank traversing a desert terrain or 
a speedboat traveling in an open sea. 

C. THESIS OVERVIEW 

For this thesis, the focus is on both scenario type, particularly 
weaponeering assessment against specific type of land vehicles moving along a 
certain stretch of road and small fast-moving boats in the open sea. An overview 
of the weapon assessment process is illustrated in Figure 1. 
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Figure 1. Overview of Weapon Assessment Process Flow 

Ideally, to perform a realistic weapon assessment, the ability to precisely 
predict the exact position of the target of specific vehicle type will enable the 
most accurate simulation of the actual combat scenario, and thereby generate 
the weapon metrics to best represent the likely outcome with the highest fidelity. 
Nevertheless, given the large amount of time required for the weapon In-Flight 
Simulation (IFS) model to produce the simulation results for each and every type 
of weapon for a particular mission, it would be impractical to have this carried out 
during the actual mission phase when time is limited. 

Therefore, for the case of the land target, it is necessary to categorize the 
wide distribution of roads into discrete road types from which a random road, 
which is representative of each of the road types, is used in the weapon IFS 
model for the full weapon assessment to be carried out. Similarly, for the case of 
the sea target in an unconstrained environment, it is necessary to predict the 
likely evasive maneuvers that the boat will adopt in its approach, and carry out 
the full weapon assessment in the weapon IFS model. This forms the first part of 
the thesis. 


3 


















Due to the vast number of roads available, it is nearly impossible to carry 
out the full weapon assessment whenever a new target emerges. As such, a 
generalization of the types of roads available is required, and a representative 
road of each type is to be generated and analyzed. In Chapter II a discussion of 
the generation of synthesized roads of different types based on type of road 
profile expected for each type is offered for the land target. The discussion 
addresses how the various roads are categorized and how the representative 
road for each category is generated. For the sea target, the route that a boat is 
likely to follow in its attempt to evade weapons during its escape or attack is 
designed and generated using an existing realistic boat simulation model. 

In Chapter III, we examine how the velocity profile of the land vehicle 
traveling along each of the synthesized roads generated can be computed based 
on the road profile. It is important to model the position and velocity profile of the 
target vehicle as closely as possible to its actual motion to ensure that the 
simulation results from the IFS model will be nearly identical to the actual 
weapon performance under actual combat scenario. For the sea target, the 
velocity profile is generated as part of the simulation. Further processing of the 
data generated is required, however, to transform the data points into the Time 
Space and Position Information (TSPI) format necessary for the weapon IFS 
model. 

Thereafter, Chapter IV contains a study of how the position of both the 
land and sea vehicles, based on the velocity profile computed In Chapter III, can 
be determined and generated at discrete time intervals, so that the TSPI can be 
fed into the weapon IFS model for the weapon simulations to be carried out. The 
IFS model requires the TSPI Inputs for the target to be at one-second Intervals, 
along with the velocity profile during each time-step. 

The weapon IFS model simulations will not be part of this thesis due to the 
classified database from which it draws the weapon information. The 
weaponeering community will make use of the TSPI Inputs to generate the 
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associated simulation results. Similarly, the simulation results are classified and 
will not be disclosed in this thesis. 

The second part of the thesis then focuses on enabling the users to make 
use of simulation results generated earlier for making weaponeering assessment 
during the actual mission planning. For this part of the thesis, the focus is on the 
land target scenario, due to the unconstrained nature of the sea target scenario. 

Chapter V presents the digitization process of an actual route of interest 
that a target is likely to travel on. It discusses how the map selection process can 
be carried out, and how the route information can be imported into a program 
and processed into a digital form that can be used for further target velocity 
profile analysis. 

Chapter VI then makes use of the digitized route to compute the likely 
velocity profile of the target, using such information to determine the route profile 
and classify the route of interest into one of the road type classification discussed 
in Chapter II. By classifying the route into one of the road types, the associated 
weapon metrics for the particular road type can then be referred to for a close 
estimate of the weapon metrics obtained from the elaborate weapon simulations 
already carried out for the specific road type. 

Chapter VII summarizes the conclusions from this study and offers 
recommendations for further studies. 

D. LITERATURE REVIEW 

In this section, we discuss the results of the review carried out on existing 
literature for the following areas to aid the study. 

1. Road Profile Derivation 

To compute a realistic velocity estimation for a vehicle traveling down a 
particular stretch of road, it is necessary to relate the maximum design speed 
limit with the road curvature along various stretches of the road. And to do that, it 

is necessary to compute the curvature at various points along the entire route. 
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Andersson and Aronsson [4] discuss the road shape modeling concept and the 
rationale for the selection of the spline interpolation of the points. The splines 
constructed are approximations, however, and the splines may not completely 
follow the exact path of interest despite the degree of shaping done. Hence, an 
alternate method of determining the road profile is required. 

With the digitization of the route, the route will be well defined by many 
points along the entire stretch of the path. Thus, an exact method of finding the 
curvature at each and every point along the path can be used. By evaluating the 
circle that passes through the point and its immediate two adjacent points, it is 
possible to determine the curvature of the road. This will be discussed further in 
Chapter III. 

Once the curvature of the path has been obtained, the speed limit at each 
point along the road can then be computed. Donnell et al. [5] discuss the 
relationship between the horizontal road curvature and the superelevation. By 
considering the centripetal acceleration of the vehicle when it travels along a 
curved road, as shown in Figure 2, the velocity of the vehicle, as a function of the 
road curvature, the side friction factor and the road superelevation, can be 
determined. 



Figure 2. Centripetal Acceleration of Vehicle Traveling along Curved 

Road. Adapted from [5]. 
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Resolving forces to both the horizontal and vertical directions, and equating 
forces in the horizontal direction only, we have 


Fcos^ +A^sin^ = 


mV^ 

R 


( 1 . 1 ) 


Similarly, for the vertical direction, 

N cos 6 -W - F = 0 . 

The relation between F and R is given by 

F = juN. 

Using Equation (1.3), and equating R in both (1.1) and (1.2), 


IV 


mV^ 


Solving for V, 


cos^-//sin^ /?(//cos^ + sin^) 




I-jU tan 0 

With Equation (1.5), the following equation is obtained: 


Zl 

157 ? 


= 0.0k + / 


( 1 . 2 ) 

(1.3) 

(1.4) 

(1.5) 

( 1 . 6 ) 


where 

e = rate of superelevation, in percent 
f = side friction (demand factor) 

V = vehicle speed, in mph 
R = radius of curve, in feet 

With Equation (1.6), the design speed, limited by the road curvature, can 
thus be determined and incorporated as one of the constraints on the vehicle in 
subsequent velocity profile analysis. 
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2. High-Speed Boat Maneuver Estimation 

In order to more accurately predict the evasive maneuvers that the boat 
will adopt in its approach, it will be essential to first determine how maneuverable 
the vessel is. The maneuverability will ultimately determine the type of turns the 
vessel is capable of, as well as the time required and associated speed 
attainable with such turns. 


Aslan [6] discusses the maneuvering characteristics of some of the fast 
craft employed in different parts of the world. Incorporating the correlation 
derived by Lewandowski [7] between the fast craft turning radius with the length 
LWL, volume V, volumetric Froude number Fy and rudder angle S, based on the 
empirical data from the works of Denny and Hubble [8], the following equation 
was developed. 


K 

LWL 


==0.5 


1.7 + 0.0222F^ 


r \ 

2.85 


LWL ' 


1 

1 V3 j 


UJ 


(1.7) 


The Froude number, in turn, is derived from the vessel approach speed 
and the volume V through the following equation. 


F,= ( 1 , 8 ) 

The normalized turning radius ^ computed and the approach speed 

can then be used to derive the speed during the turn 1]^ through the following 
equation by Denny and Hubble [8]. 


U=U,. 


l + K^ 


f \2 

‘ LWL^ 


V y 


(1.9) 


The empirical constant can be obtained using the Froude number Fy 
and rudder angle S using the following equation. 
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2 


( 1 . 10 ) 



To ensure the stability of the vessel during a turn, the Nomoto K-T model 
is studied in [6]. To simplify the model, the first order transfer function by Nomoto 
was adopted. When expressed in time domain, the following equation is derived 
and used in the model. 


where 


Tr + r = KS. 


r = y/. 


( 1 . 11 ) 

( 1 . 12 ) 


Assuming the yaw inertia coefficient T to be zero, Aslan [6] then derives 
the following relation between the turning ability K and the vessel approach 
speed Ua, the turning radius Rc and the max rudder angle S^^ax- 




U„ 


R.xS^ 


(1.13) 


Using the computed value of K, substituting Equation (1.13) into Equation 
(1.11) and introducing a correction factor rj of value 0.909 to ensure consistency 
with the turning radius Rc in the Lewandowski’s equation [7], the following 
equation is obtained. 

i^ = KSij. (T14) 


Equation (1.14) can then be integrated with respect to time for the 
computation of the velocity profile using the following equations. 


x = U^ sin y/ 
y = U^ cos y/ 


(1.15) 


Further integration of Equation (1.15) can generate the position profile of 
the vessel with respect to time. 
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3. Map Coordinates Transformation 

To digitize a road, it is necessary to convert map locations in the Geodetic 
Coordinates System to that of a local reference frame. Both Farrell [9] and 
Oriels [3] discuss the conversion of the Global Positioning System (GPS) 
readings from the WGS-84 geodetic system to a local East-North-Up (ENU) 
Cartesian axis. 

The basic reference frame for terrestrial navigation is the Earth-Centered- 
Earth-Fixed (ECEF) geodetic coordinate system, where the coordinates are in 
latitude 0, longitude X, and altitude h. This is also referred to as the global 
reference frame, which is most commonly associated with the GPS. 

An alternate reference frame is known as the ECEF rectangular system, 
where the x-axis passes from the origin at the Earth’s center to the prime 
meridian, the z-axis passes through the North Pole, and the y-axis orthogonal to 
these two axes. This frame is convenient, as it allows for easy relation to the 
ENU frame that is used subsequently. 

To use the ECEF geodetic system, the shape of the Earth’s ellipsoid 
needs to be defined, and the parameters for the assumed shape used in GPS, 
known as the World Geodetic System - 1984 (WGS-84), are defined by the 
equatorial and polar radii, a and b, of 6,3778,137.0 m and 6,356752.3142 m, 
respectively. The ellipse is presented in Figure 3. 
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Figure 3. WGS-84 Ellipse 


The ellipsoid flatness is defined by 

f = "—. (1.16) 

a 

The eccentricity of the ellipsoid is defined by 

e = ^f(2-f). (1.17) 

The distance from the surface of the ellipsoid to the z-axis along the 
ellipsoid normal, as a function of 0, is defined by 

I ■ ( 1 . 18 ) 

yjl-e sin (j) 

The ECEF rectangular coordinates can then be calculated from the 
geodetic system using the following: 


X = (/z + A^) cos (^cos A., 

(1.19) 

y = (/z + A^) cos sin A , 

(1.20) 


11 













( 1 . 21 ) 


Thereafter, in order to form a right-handed coordinate system that is 
analogous to the Cartesian coordinate system, it is essential to further transform 
from the ECEF rectangular system to that of the END system. This rectangular 
system is defined with the origin at a point of interest on the Earth’s surface, with 
the vertical axis along the ellipsoid normal, and the other two axes orthogonal 
and aligned to the local geographic East and North directions. The END is 
presented in Figure 4. 



J'ecef 


Figure 4. ECEF Geodetic and ENU Coordinate Frame. Source: [3]. 

With the origin of the ENU defined as (xo,yo-^o) in the ECEF rectangular 
frame, to transform the ECEF rectangular frame to the ENU frame would require 
a translation of the origin to (xo,yo-^o). followed by a rotation of the axes to align 
with those of the ENU. The complete transformation is given by 
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( 1 . 22 ) 




-sin/I 

cos/I 

0 

X-Xg 

n 

= 

- sin ^ cos A 

- sin ^ sin A 

cos^ 

o 
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cos, (j) cos X 

cos ^ sin A 

sin^ 

_Z-Zo 


With the transformation in Equation (1.22) into the ENU system, which is 
analogous to the familiar Cartesian system, the digitized route information can 
then be easily further processed in subsequent analysis. 
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II. ROUTE GENERATION 


A. LAND ROUTE GENERATION 

In this section, the aim is to generate routes that are representative of the 
different road types. Before we begin to generate synthesized road for each road 
type, it is essential to first define a quantitative approach to classify the road 
types. Thereafter, we will generate the roads based on the characteristic 
parameter associated with each road type. 

1. Statistical Road Type Classification 

Driels [10] introduced the concept of road type characterization through 
the analysis of actual roads. A statistical analysis is carried out to study both the 
heading and change in heading for four distinct types of real roads to determine 
their respective means and standard deviation. The statistical results, in 
histogram form, are presented in Figures 5 to 8. 
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Figure 5. Statistical Analysis of Coast Ridge Road. Source: [10]. 
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Figure 6. Statistical Analysis of Carmel Valley Road. Source: [10]. 
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Figure 7. Statistical Analysis of Flighway 1. Source: [10]. 


16 

























Hwy 101 


Heading Mean (deg) 138.583 
Heading StDev (deg) 22.086 
A Heading Mean (deg) -0.009 
A Heading StDev (deg) 0.434 
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Figure 8. Statistical Analysis of Highway 101. Source: [10]. 

The results of the preceding statistical analyses are consolidated in Table 
1 for ease of comparison. 


Table 1. Statistical Analyses of Actual Roads 


Road type 

Mean 

heading 

Sigma 

heading 

Mean 
change in 
heading 

Sigma 
change in 
heading 

Coast Ridge 
Road 

121 

54.6 

-0.02 

13.42 

Carmel Valley 
Road 

120 

43.3 

-0.01 

4.93 

Pacific Coast 
Hwy (Hwy 1) 

138 

31.0 

-0.10 

3.63 

Highway 101 

139 

21.1 

-0.01 

0.49 
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From the results obtained, it can be easily observed that the change in 
road heading generally follows a normal distribution, and the standard deviation 
in the change of heading increases with the number of curves and degree of 
curvatures expected on the roads. From this, for the synthesized road 
generation, the synthesized roads are assumed to follow a normal distribution 
and parameterized by the standard deviation of the road change in heading. A 
total of six different road types are identified, and their associated standard 
deviations in change of heading are summarized in Table 2. 


Table 2. Road Type Classification by Change in Road Fleading 


Road 

type 

Sigma change in 
heading 

Mean sigma change in 
heading 

1 

0< to <1 

0.5 

2 

1< to <4.5 

2.5 

3 

4.5< to <9.5 

7 

4 

9.5< to <14.5 

12 

5 

14.5< to <19.5 

17 

6 

19.5< to <30.5 

25 


2. Route Synthesis by Road Type 

Using the parameters defined in Table 2 alone is not enough to fully 
characterize and produce synthesized roads that are representative of the 
original road. A comparison of the original Carmel Valley and the synthesized 
road reproduced using the same road parameters is presented in Figure 9. 


18 






Figure 9. Comparison of Carmel Valley Road and 
Synthesized Reproduction 

It can be clearly seen that the synthesized reproduction of the Carmel 
Valley road is totally random and lacks a general direction of travel. Therefore, 
synthesizing roads in such an unconstrained manner is not suitable, and 
additional constraints are required to produce synthesized roads that can be 
more representative of the original. To prevent excessive randomness in the path 
generation and to impose a general road heading, the mean heading and the 
standard deviation in the heading are included for the road synthesis. 

The road synthesis is carried out by iteratively generating small segments 
of path lengths with change of heading randomly drawn from a normal 
distribution with the mean and parameters as defined in Table 3. These road 
segments are then incrementally added on to an arbitrarily chosen starting point 
until a finite route distance is achieved. In addition to this, a further constraint is 
also imposed on the road heading by restricting the maximum change in heading 
to no more than two standard deviations of the heading. This is illustrated in 
Figure 10. 
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Table 3. Final Road Type Classification 


Road 

type 

Mean heading 

Sigma 

heading 

Sigma change in 
heading 

Mean sigma change 
in heading 

1 

135 

25 

0< to <1 

0.5 

2 

135 

29 

1< to <4.5 

2.5 

3 

135 

41 

4.5< to <9.5 

7 

4 

135 

53 

9.5< to <14.5 

12 

5 

135 

66 

14.5< to <19.5 

17 

6 

135 

73 

19.5< to <30.5 

25 



Figure 10. Route Synthesis Path Constraint 


With these additional constraints imbued into the program, a new 
synthesized road is generated for verification of the updated program. The new 
synthesized road, while not exactly the same as the actual Carmel Valley road, 
which is expected due to the randomness in the generation of paths during the 
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road synthesis, is much closer to the actual road than before, and more 
representative of a road in the specific road type. The roads are provided in 
Figure 11 for comparison. 



Figure 11. Comparison of Digitized Carmel Valley Road with New 

Synthesized Road 


With the new changes incorporated, and having verified the ability of the 
program to produce a synthesized road that is reasonably representative of 
actual digitized road, the representative synthesized roads for each of the six 
distinct road types are generated and presented in Figures 12 to 17. The 
Cartesian coordinates of all the points making up the newly synthesized roads 
are saved for subsequent work in Chapter III for conversion into the TSPI input 
format required for weapon performance analysis by the IFS model. 

As intended, the roads generated all have a common heading towards the 
southwestern direction, and the road profile gets increasingly winding as the 
standard deviation in the change of heading increases. From the plots of the 
routes generated, it can be easily observed that the routes are as per design, 
and they are ready to be worked on further in subsequent chapters. A copy of the 
MATLAB codes for the generation of routes for the various road types can be 
found in Appendix A. 
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Figure 12. Synthesized Road for Type 1 


Type 2 



Figure 13. Synthesized Road for Type 2 
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Figure 14. Synthesized Road for Type 3 


Type 4 



Figure 15. Synthesized Road for Type 4 
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Easting (km) 

Figure 16. Synthesized Road for Type 5 


Type 6 



Easting (km) 

Figure 17. Synthesized Road for Type 6 
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B. SEA ROUTE GENERATION 


In earlier discussion, the target’s route of travel is constrained by a specific 
path, and these paths can in turn be characterized by their profile into different 
road types, from which a representative path can be generated for further 
analysis. Yet, for a vessel moving in the open sea, there is no such constraint to 
restrict it within a fixed path. Instead, it is free to move wherever it wants, 
whenever it wants. In order to generate a realistic prediction of the target’s path 
of travel on the open sea, certain assumptions about its travel are required. 

For this study, it is assumed that the general heading for the target, when 
viewed on an X-Y Cartesian plane, is towards the right, along the x-axis. To 
reach its target in the shortest amount of time, it always maintains its approach 
speed at the maximum possible. To try to evade a weapon directed at it, though, 
the vessel must execute turns at the highest turning speed possible, while 
maintaining its general motion towards the desired heading. This is termed as the 
desired travel path for the sea vessel, as illustrated in Figure 18. 



Figure 18. Sea Vessel Predicted Path of Travel with Evasive Maneuvers 
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Of course, a vessel cannot perform sharp turns, as depicted by the 
desired travel path in Figure 18, due to the turning radius required when making 
turns at high speeds. As such, a more realistic path for a sea vessel executing 
evasive maneuvers at high speeds is more accurately represented by the blue 
path, as shown in Figure 18. 

Furthermore, for the turning angle 6, the sharper the angle, the higher is 
the speed of the sea vessel along the x-direction. Thus, a large angle will enable 
the target to move down the x-axis in the shortest possible time while executing 
the evasive maneuvers. A large angle turn, however, will not be aggressive 
enough for an evasive maneuver. Flence a sharp turn angle of 60 degrees is 
used in the predicted route of the sea vessel to simulate more aggressive 
evasive maneuver while maintaining a substantial speed along the x-axis for the 
purpose of this study. 
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III. VELOCITY PROFILE COMPUTATION 


A. VELOCITY PROFILE FOR LAND VEHICLE 

In this section, the aim is to determine the velocity profile of the land 
vehicle traveling along a constrained path. The model maximum speed of the 
vehicle is constrained by both the vehicle maximum speed and the road design 
speed. As mentioned in Chapter I, in order to make use of Equation (1.6) to 
compute the maximum road design speed limit along a curve, it is necessary to 
determine the radius of curvature. 

1. Radius of Curvature Computation 

Six distinct routes, made up of thousands of intermediate points with their 
associated Cartesian coordinates, have been generated earlier. Due to the large 
number of points in each route, it is assessed that it will be computationally less 
intensive to solve for the radius of curvature at a point using a circle that passes 
through the point and the two points immediately adjacent to it, as shown in 
Figure 19. Furthermore, this method produces an exact solution for the radius of 
curvature, while there may be certain loss of accuracy, even with good shaping, 
with the method of spline interpolation, as suggested by Andersson [4]. 


(X2,y2) 



Figure 19. Radius of Curvature Computation Using Three Known Points 
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The gradients of the lines connecting the first and second pair of points 
are given by 


= 


3^2-Ji 


-X 


^3-^2 


X, -x. 


(3.1) 


The equation of the lines connecting the midpoint of line of gradient rtia 
and that of gradient mi, to the center of the circle of curvature are given by 




Jl + ^2 __1 
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Equating the value of using the equation of lines in (3.2), we get 
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(3.3) 


Using Equations (3.2) and (3.3), the center of the circle of curvature, and 
hence the radius of curvature, can be computed for every point on the 
synthesized route, less the start and end points. Having determined the radius of 
curvature, we can then commence on establishing the velocity profile of a target 
traveling along the designated route. 


2. Maximum Model Speed Limit Computation 

Having computed the radius of curvature in the previous section, we can 
then proceed to compute the maximum road design speed at each point of the 
route based on the road profile using the Equation (1.6) discussed by Donnell [5], 
as mentioned in Chapter II. Ideally, the vehicle should travel at the road design 
speed limit in order to minimize the time taken to reach the destination. However, 
this is only possible if either the route that the vehicle is traveling on is fairly 
straight with little change on the maximum speed limit, or the vehicle has 
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extraordinary acceleration and braking capabilities to enable it to change to the 
desired speed at will. Even a Formula 1 race driver racing down a straight stretch 
of road must brake hard and early to successfully negotiate a sharp bend ahead 
without skidding off the road, and thereafter accelerate hard down the straight 
road to reach maximum speed. Therefore, to calculate the realistic maximum 
model speed that the vehicle is able to travel along any route, not only is it 
important to determine the road design speed limit to prevent the vehicle from 
skidding out during a turn, it is also equally important to take into consideration 
the vehicle’s performance capabilities, i.e., the vehicle’s maximum speed, as well 
as its acceleration and deceleration limits. 

To determine the maximum model speed for the vehicle at any point along 
the route of travel, it is necessary to compare the design road speed limit at each 
point with the maximum vehicle speed limit, and choose the lower of the two. 
This will set the maximum model speed that the vehicle cannot exceed at every 
point along the road. 

For the purpose of weaponeering assessment, there are generally three 
classes of vehicles that are of interest. For this thesis, the focus is on the first 
class of vehicle where the maximum vehicle speed is expected to be around 60 
mph. For the vehicle acceleration limit, the time taken for an average vehicle to 
reach 100 km/h is around ten seconds, and assuming a linear acceleration 
profile, the acceleration limit is assumed as 0.3g for our study. Assuming that the 
braking is twice as effective, the deceleration is assumed as 0.6g. The 
breakdown of the expected vehicle performance for all three classes is 
summarized in Table 4. With the speed, acceleration, and deceleration limits all 
established, it is now possible to compute the actual speed of the vehicle at 
every point along the route of interest. 
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Table 4. Types of Land Vehicles of Interest 


Vehicle type 

Max speed 
ms'^ (mph) 

Max 

acceleration (g) 

Max braking (g) 

Fast 

26.82 (60) 

0.3 

0.6 

Slow 

13.41 (30) 

0.15 

0.3 

X-slow 

6.71 (15) 

0.08 

0.16 


To determine the actual speed of the vehicle at each point along the route, 
it is important to first identify the points near the current point where the model 
maximum speed is lower than that for the current point. Imagining ourselves as 
the driver of the vehicle and assuming we reach the start of the route of interest 
with a certain reasonable starting speed, at the current instance, we are only 
concerned with any point lying ahead that requires us to take on a lower speed 
than our current speed, as these points will determine the time and amount of 
braking needed to reach the lower speed required in time. 

An illustration of the deceleration constraint is provided in Figure 20. The 
vehicle is traveling at 27 m/s currently, and a point 60 meters ahead has a model 
maximum speed of 12.4 m/s, and the vehicle is capable of decelerating at 6 m/s^. 
Assuming a linear deceleration profile, the relationship between the distance 
traveled s, the initial velocity u, and the final velocity v is given by 

v^=u^ + 2as. (3.4) 

Using Equation (3.4), the vehicle is capable of reaching 12.4 m/s from 27 m/s in 
a distance of 48 meters. Thus, there is no need for the vehicle to reduce speed at 
the current instance, and it can continue traveling at the model maximum speed 
until point 2. If the model maximum speed at point 2 is higher, it can even 
accelerate to 27.4 m/s before decelerating to 12.4 m/s at point 7. Once it reaches 
point 2, the vehicle can maintain the model maximum speed, and it will need to 
begin decelerating in order to meet the model maximum speed at point 7 


30 








subsequently. For point 3, due to the lower model maximum speed at point 7, the 
maximum speed possible is only 25.2 m/s, and the vehicle will need to 
decelerate to this speed when it reaches point 3. 



Figure 20. Illustration of Deceleration Constraint in Speed Computation 

Thus far, the focus of the example using Figure 20 has been on the effect 
of the lowest model maximum speed of 12.4 m/s at point 7 on the maximum 
possible speed at all the points before point 7. Nevertheless, it is important to 
note that the effects of this for all the points with lower model maximum speeds 
will need to be considered instead. A model maximum speed of 12.5 m/s at point 
6, though still higher than that at point 7, would have a greater impact on the 
traveling speed at the earlier points due to the shorter stopping distance 
available. Thus, the model, as it steps through each point on the path, will need 
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to consider all the points with lower model maximum speeds within a specific 
lookahead distance, and compute the maximum possible speed for the vehicle at 
the current point based on the deceleration limit and the stopping distance 
available to each of these possible deceleration “limiting factors.” 

Along the same line of reasoning, the application for the acceleration limit 
constraint must also be considered as the vehicle travels from point to point 
sequentially along the route of interest. For the acceleration constraint, the points 
of interest are those prior to the current point with lower model maximum speeds. 
An illustration of the acceleration constraint is provided in Figure 21. which is a 
continuation of the previous example in Figure 20. 



Figure 21. Illustration of Acceleration Constraint in Speed Computation 
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For point 8, the maximum speed possible due to the acceleration limit is 
14.6 m/s; thus, the constraint in this case is due to the model maximum speed. 
For point 9, the maximum speed of the vehicle at point 7 is only 12.4 m/s, and 
the vehicle is capable of accelerating at 3 m/s^. Assuming a linear acceleration 
profile, using Equation (3.4), the maximum possible speed due to the 
acceleration limit from point 7 over the 20m distance is 16.5 m/s, which is lower 
than the model maximum speed of 17.6 m/s. Yet, the realizable speed of the 
vehicle is lower at 16.0 m/s instead due to the effect of the model maximum 
speed at point 8, with a pick-up distance of 10 m only. 

Thus, the model, as it steps through each point on the path, will need to 
consider all the points with lower model maximum speeds within a specific 
lookback distance and compute the maximum possible speed for the vehicle at 
the current point based on the acceleration limit and the pick-up distance 
available to each of these possible acceleration “limiting factors.” 

Flaving introduced both the deceleration and acceleration constraints for 
the model, we then need to combine both these constraints for the computation 
of the actual speed of the vehicle. This is easily achieved by comparing the 
maximum speeds possible due to both the acceleration and deceleration 
constraints with the model maximum speed at every point along the path, and 
choosing the lowest of the three to satisfy all the constraints. All three constraints 
are presented for the two earlier examples in Figures 22 and 23. 

From Figure 22, we can see that the actual speed of the vehicle at points 
1-8 is primarily constrained by the deceleration limit and the model maximum 
speed. This is expected as the model speed limit along this stretch is decreasing 
in general. 
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Figure 22. Computation of Speed under All Three Constraints 

for Points 1-8 



Figure 23. Computation of Speed under All Three Constraints 

for Points 7-14 


34 



















From Figure 23, we can clearly see that the speed for the vehicle at points 
9 to 11 is determined by the acceleration limits of the vehicle. The computation of 
the speeds at points 10 and 12 is of interest, as the values of the speed due to 
the acceleration, deceleration, and model maximum speeds constraints are all 
distinctly different. To meet the constraints of all three limits, the speed of the 
vehicle is chosen to be the lowest of the three limits. For point 10, the deciding 
“limiting factor” is that of the acceleration limit, while that for point 12 is the 
deceleration limit determined by the low model speed limit at point 13. 

The approach that we have discussed thus far is then applied to all the 
points along the path of interest to generate the highest possible speed of the 
vehicle at each of the points, thereby generating the velocity profile of the vehicle 
that satisfies all three physical constraints discussed earlier. 

The velocity and acceleration profile for each of the six road types are 
plotted along the entire path to verify that the speed and acceleration constraints 
are fully satisfied. Due to the large number of points involved for each plot, a plot 
of 100 sampling points is provided for each of the road types for better clarity. 
The respective plots are provided in Figures 24 to 35. 
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Figure 24. Speed Profile for Road Type 1 




Figure 25. Acceleration Profile for Road Type 1 































Figure 26. Speed Profile for Road Type 2 



Figure 27. Acceleration Profile for Road Type 2 
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Figure 28. Speed Profile for Road Type 3 



Figure 29. Acceleration Profile for Road Type 3 


38 








































































































Figure 30. Speed Profile for Road Type 4 



Figure 31. Acceleration Profile for Road Type 4 
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Figure 32. Speed Profile for Road Type 5 



Figure 33. Acceleration Profile for Road Type 5 
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Figure 34. Speed Profile for Road Type 6 



Figure 35. Acceleration Profile for Road Type 6 
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From all the velocity and acceleration plots presented, we can see that the 
speeds of the vehicle are kept at less than or equal to the model maximum 
speeds, while the acceleration of the vehicle does not exceed the limits imposed. 
This implies that the solution meets all the kinematic constraints satisfactorily; as 
such, the results obtained from the approach discussed earlier are valid. 

During the course of the thesis, an optimization approach is explored to 
determine the solution for a time-minimization problem. Gill [11] discusses the 
use of the SNORT software for the constrained optimization problems, which is 
applicable to this thesis. Numerous optimization runs are attempted with no 
meaningful results obtained until the end of the thesis work. As further tuning of 
the optimization program is required to generate the results that meet all the 
constraints, the optimization approach is brought to a halt, and further exploration 
of this area is recommended for future work to obtain the optimal solution to the 
time-minimization problem. 

An alternate approach for the velocity computation is also explored toward 
the end of the thesis to include the effects of the vehicle dynamics in the 
computation of the vehicle velocity profile. Kaminer [12] developed a model in the 
Simulink environment to include the effects of vehicle dynamics on the vehicle 
trajectory and its velocity profile as it travels down a particular route. 

The velocity profile computation is carried out using the model developed 
by Kaminer [12], and the results for the initial approach and that obtained using 
the Simulink model are provided in Figure 36. Due to the large number of data 
points involved, a small sample of the data points is illustrated in the plots for 
ease of comparison. 
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Figure 36. Comparison of Velocity Profile 


From the obtained, it can be seen that the velocity profile generated by 
Simulink does not produce the same sharp changes in velocity that are observed 
in the velocity profile obtained initially; as such, the Simulink model is able to 
better capture the effects of vehicle dynamics and produce a more realistic 
representation of the actual vehicle velocity profile. However, it can also be 
observed that the velocity profiles generated in both approaches are largely 
similar. Therefore, the velocity profile generated using the approach discussed 
earlier provides a very good approximation of the actual vehicle performance 
under the given constraints. 

In view of the limited time remaining for the completion of the thesis work, 
the Simulink model approach is not studied further. As the results from the initial 
approach are very close to those obtained from the Simulink model, the original 
set of velocity profiles computed is retained. For the rest of the thesis, the focus 
is on using the results from the initial approach used for the computation of the 
vehicle velocity profile. For future work in similar areas, the Simulink model 
approach is strongly recommended. A copy of the Simulink model for velocity 
computation is provided in Appendix B. 
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B. VELOCITY PROFILE FOR A SEA VESSEL 


For the case of the sea vessel, the approach is different, as the simulation 
model by Aslan [6] is able to generate the velocity profile of the sea vessel 
throughout the course of its motion. The course of the sea vessel is determined 
by the rudder inputs, and its maneuverability dependent on the sea vessel 
characteristics. The sea vessel characteristics that are required as inputs to the 
simulation model for two examples of sea vessels of interest are summarized in 
Table 5. Therefore, the aim is to generate a suitable set of rudder input to 
produce the desired path of travel, shown In Figure 18, for the sea vessel and 
generate the associated velocity profile. 


Table 5. Inputs Characteristics of Sea Vessels. Adapted from [13], [14]. 


Sea Vessel 

Length on 
Waterline 
(m) 

Displacement 

(ton) 

Maximum 

Rudder 

Angle 

(deg) 

Approach 

Speed 

(knots) 

Bladerunner 51 

15 

16 

30 

65 

Boston Whaler 

7.75 

5 

30 

45 


To generate the correct rudder inputs to produce the desired path of travel 
for the sea vessel, it is necessary to first understand the maneuverability of the 
sea vessel, which is vessel specific and dependent on the vessel length, weight, 
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and maximum speed. To do that, we first set the rudder angle input at the 
maximum rudder angle throughout the simulation run. This will simulate the sea 
vessel making the tightest turn when going at its maximum speed, and from this, 
we can determine the time taken for the boat to make a complete turn. The 
results of the maneuverability study are shown in Figure 37. 


Sea Vessel Maneuverability Study 



Figure 37. Sea Vessel Maneuverability Study 

Through the maneuverability study, we determine that the Bladerunner 51 
takes around 41.05 seconds to turn through a complete round of tactical 
diameter 389 meters at its maximum turning speed of 29.8 m/s. By contrast, the 
Boston Whaler takes around 13.04 seconds to turn through a complete round of 
tactical diameter 68 meters at its maximum turning speed of 16.3 m/s. With this 
information, we can then proceed to work on generating the rudder inputs for 
producing the path of travel desired. 

Since the course of the sea vessel is driven by the rudder inputs, it is 
necessary to compute the duration of the rudder inputs as the sea vessel goes 
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through each turn. This can be easily done by considering the angle that the ship 
turns through, since we have already determined the time it takes for the ship to 
turn through a complete round. Referring to the angle 9 in Figure 18, the angle 
that the sea vessel turns through can be computed by computing 180°-6>, 
based on the geometrical relationship shown in Figure 38. 



Figure 38. Geometrical Relationship for Sea Vessel Turn Angle 

The turn angle computation allows us to compute the turn duration for 
every turn in the path. Putting all the turn durations into a combined rudder input 
for the simulation, the control to the system is presented in Figure 39, and the 
resulting output corresponding to the input is shown in Figure 40. 
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Rudder Angle 



Figure 39. Sea Vessel Rudder Angle Inputs 


Course of Sea Vessel on X-Y Plane 



Figure 40. Sea Vessel Simulated Course of Travel 
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From the results obtained, it can be clearly seen that the sea vessel 
responds as per intended, turning with the changes in the rudder angles, moving 
the sea vessel along the sea route that was planned and discussed in the 
previous chapter. 

Having achieved the desired path of travel for the sea vessel, and with the 
associated velocity profile generated by the simulation model, we are now ready 
to convert the position and velocity profile into the TSPI format required for the 
weapon IFS model. A copy of the Simulink model for the sea route simulation 
can be found in Appendix C. 
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IV. TIME DISCRETIZATION 


A. TRAVELING TIME COMPUTATION 

Having computed the position and velocity profiles of both the land and 
sea vehicle at every point on their paths of travel, we can then compute the time 
duration for the vehicle to travel between each consecutive pair of points on the 
path. This can be done using the kinematic model with linear acceleration, as 
shown in Figure 41. 



Figure 41. Kinematic Model for Time Duration Computation 

With the distance between each pair of consecutive points known, the 
time duration can be computed using the following equation. 


At 


i — 




(4.1) 
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After we have computed the time duration the vehicle spent traveling 
between points, we can then proceed to recreate the path of travel with points at 
discrete time intervals to generate the data for the TSPI inputs. 

B. CONVERSION TO FIXED-TIME INTERVAL 

To recreate the path of travel with points at discrete time interval, the 
approach is to iteratively step through each point in the path, sum up the time 
duration for the vehicle to travel to the current point, and check whether the sum 
exceeds multiples of one second. If the total time duration is less than one 
second, the model will proceed with the next iterative step to the next point on 
the path. If the total time duration is more than one second, a new point on the 
new path with a discrete time interval is created by interpolating between the 
current point and the previous point for which the time duration is exactly one 
second, and then subtracting one second from the total time duration. This 
continues until the total time duration is less than one second, and the model will 
proceed with the next iterative step to the next point on the path. A visual 
representation of the iterative process is presented in Figure 42. 



Figure 42. Creation of New Points at Discrete Time Intervals 
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Besides creating the new points at every discrete time interval, it is 
important to compute the velocity of the vehicle at each of the newly created 
points as well. This can be done during the same iterative process by 
interpolating the velocity of the new point from the velocities of the vehicle at the 
points between which the new point is being created, based on the position of the 
new point relative to the positions of the points for the interpolation. This 
approach is acceptable due to the linear acceleration model that has been 
adopted for the model earlier. The general equation for the interpolation of the 
velocity of the new point j from those of point i and (i + 1) is given in 
Equation (4.2). 


U I = M; + 

j * 


At.^ 




(4.2) 


Having derived the position and velocity profile of the new points for the 
path at each discrete time interval, the required inputs for the TSPI file can now 
be generated. 


C. NEW DATA GENERATION FOR LAND ROUTES 

Using the approach that has been discussed, the routes that we have 
generated previously can now be converted to ones with the discrete time 
intervals. The new paths for the routes of road type 1 to road type 6 are created 
and presented in Figures 43 to 48, respectively. 
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Figure 43. Type 1 Route at Discrete Time Intervals 
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Figure 44. Type 2 Route at Discrete Time Intervals 
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Figure 45. Type 3 Route at Discrete Time Intervals 
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Figure 46. Type 4 Route at Discrete Time Intervals 
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Figure 47. Type 5 Route at Discrete Time Intervals 


Type 6 Road 



Figure 48. Type 6 Route at Discrete Time Intervals 
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Due to the large number of data points used for the generation of the six 
routes for type 1 to type 6, only 5 percent of the points are presented in order to 
see the new data points on the plots more visibly. Plotting the entire range of 
data points will only generate plots that look similar to those produced initially. 

From the plots obtained, it can be observed that the points get closer 
together along a curvy portion of the route, and further apart when the route gets 
straighter. As the vehicle travels faster along a straight road, the amount of 
distance covered with each second of time is greater. Similarly, when the road 
gets curvier, the vehicle speed is reduced, and the amount of distance that can 
be covered in a second of time is reduced accordingly. As such, the results 
obtained are as expected; therefore, the results are deemed valid and 
acceptable. 

A copy of the MATLAB codes for the computation of the velocity profile 
and the conversion of the vehicle position along the route to discrete time 
intervals can be found in Appendix D. 

D. NEW DATA GENERATION FOR SEA ROUTES 

As the number of points for the path of the sea vessel is significantly 
lower, the plots of the original path and those of the discrete time interval are 
presented in Figures 49 and 50, respectively, for ease of comparison. 
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Figure 49. Original Sea Vessel Route 


Course of Sea Vessel on X-Y Plane - 1s interval 



Figure 50. Sea Vessel Route at Discrete Time Interval 
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Similar to the results obtained for the land vehicle, the data points get 
closer together as the curvature of the path gets higher, and further apart as the 
curvature of the path decreases. This is again in accordance with the 
expectations, and hence the results are valid and acceptable. A copy of the 
MATLAB codes for the conversion of the sea routes to discrete time intervals can 
be found in Appendix E. 
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V. ROUTE DIGITIZATION 


The past few chapters, which constitute the first part of the thesis, focused 
on generating the routes of interest for weapon assessment study and converting 
them into the TSPI inputs that can be used in the weapon IPS model for analysis. 
For the rest of the chapters, which constitute the second part of the thesis, the 
aim is to come up with a methodology to convert an actual route selected from a 
map and convert it into a digital form for route characterization and weapon 
assessment. An additional requirement for this methodology is a Visual Basic for 
Applications (VBA) environment. This environment enables intended users to 
utilize the program in a restricted work area with no access to advanced software 
for performing such analysis. 

A. SELECTING THE ROUTE OF INTEREST 

For the route selection, we want to take an approach that is simple, user- 
friendly, and easily accessible. With that in mind, we start off with the selection of 
the map route using Google Maps [15], a web-based electronic map, commonly 
used by many around the world. Google Maps is easily accessed using web 
browsers and an active internet connection. The limitation to this approach is that 
the routes available for selection are limited to areas that Google has managed 
to access and map out. 

Selection of a route using Google Maps [15] is simple and straightforward, 
as this can be done using the “Directions” function that can be accessed with just 
a click of a mouse button. Thereafter, the selection of the route of interest can be 
done by clicking on the locations on the map to specify the start and end points 
of the route. Google Maps [15] will then process and generate the possible 
routes to get from the start point to the end point. Then, users can choose from 
the routes that Google Maps [15] has suggested or can specify the route that 
they want by including additional waypoints. An illustration of this for a route 
along the Carmel Valley road is presented in Figure 51. 
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After selecting the route of interest, we can then proceed to extract the 
route information and convert it into a digital form that allows for further 
processing of the information that we are keen to sieve out. 

B. IMPORTING THE ROUTE OF INTEREST 

One of the common file formats for presenting map information is the GPX 
file format that GPS devices utilize. For the purpose of this thesis, we use the 
GPX file format for the extraction of the route information. While Google 
Maps [15] provides a user-friendly interface for route selection, it does not allow 
for the extraction of the route information. As such, we need to utilize other 
avenues to get the data we need from Google Maps. A quick search on the web 
will reveal many tools, available for purchase or at no cost, for performing this 
task. With the key consideration of accessibility in mind, a web-based conversion 
tool is preferred. Therefore, for this study, a third-party website, MapsToGPX 
[16], is chosen not only for its accessibility online but its simplicity of use. 
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Extraction of the map information can be easily done by inserting the 
Uniform Resource Locator (URL), also commonly known as the web address, of 
the route that we have selected in Google Maps [15] earlier and clicking the 
“Let’s Go” button to process the route information into a downloadable GPX file. 
An illustration of the usage of the map conversion tool is provided in Figure 52. 
With this file, the information of all the points that define the route selected in 
Google Maps is retrieved in digital format. 



Figure 52. Generation of Route Data in GPX Format 

After obtaining the route information in GPX file format, we must import 
the information into a software tool so that the Information can be read and 
processed easily. For this purpose, both MATLAB and Microsoft EXCEL are 
capable of importing the GPX file and presenting the map data in terms of their 
latitude and longitude coordinates in array form. With this, we are ready to 
process the information for further analysis of the route. 
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VI. ROUTE PROCESSING AND ANALYSIS 


Using the approach discussed earlier under “Map Coordinates 
Transformation” in Chapter I, we can convert the map data extracted from the 
GPX file from the geodetic latitude and longitude coordinates to the ENU 
coordinate system. In this way, the Easting and Northing coordinates are 
equivalent to the familiar Cartesian coordinate system. 

The key parameter that we want to extract from the route is the standard 
deviation of the change in heading of the route. This parameter enables us to 
classify the route as one of the six road types presented in Table 2 in Chapter II. 
To compute the standard deviation of the change in heading, first we need to 
standardize the route make-up for a fair statistical computation of the heading 
and change in heading. This is done by converting the route from one comprising 
irregular route points to one with regularly spaced route points, and for this 
thesis, the route points are spaced ten meters apart. 

A. CONVERSION TO REGULAR INTERVALS 

For the conversion of the route points to regular intervals, the approach is 
similar to that used for conversion to the fixed-time interval that we discussed in 
Chapter IV. 

To recreate the path of travel with points at equal intervals of ten meters 
apart, the process required is to iteratively step through each point in the path, 
sum up the distance traveled for the vehicle to reach the current point, and check 
whether the sum exceeds multiples of ten meters. If the total distance traveled is 
less than ten meters, the model will proceed with the next iterative step to the 
next point on the path. If the total time duration is more than ten meters, a new 
point on the new path with regularly spaced points is created by interpolating for 
the position between the current point and the previous point for which the 
distance traveled is exactly ten meters, and then subtracting ten meters from the 
total distance traveled. This continues until the total distance travelled is less 
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than ten meters, and the model will proceed with the next iterative step to the 
next point on the path. A visual representation of the iterative process is 
presented in Figure 53. 



Figure 53. Creation of New Points at Regular Distance Interval 


The general equation for the interpolation of the x and y coordinates of the 
new point j from those of point i and (i + 1) is given as 


Xj=X,+ 



( 6 . 1 ) 


y,' = y,- + 




( 6 . 2 ) 


This approach is used to convert the Carmel Valley route extracted from 
Google Maps [15] into a route made up of points spaced ten meters apart. The 
original and newly converted routes are presented in Figures 53 and 54, 
respectively. 
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Figure 54. Original Route Points Imported 



Figure 55. New Regularly Spaced Route Points 
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From the comparison of the visual plots of the two routes, we can see that 
the original route has been successfully converted into one with the route points 
regularly spaced at ten meters apart. Having established the methodology for 
standardizing the route in terms of their make-up, we can then proceed to 
conduct the statistical analysis for the converted route. 

B. CHANGE OF HEADING COMPUTATION 

To compute the heading, we need to establish the heading of each of the 
points on the path. We do this by calculating the slope of the line joining the 
points immediately before and after the point of interest. After the standardization 
of the route points earlier, the regularly-spaced route points allow this approach 
to give a good approximation for the slope of the tangent to the route at the 
current point. This methodology is applied for the computation of the slope at 
every point along the route. The computation of the slope of the tangent to the 
route and the heading at a route point is presented in Figure 56. 
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Figure 56. Computation of Heading 

With this approach, the heading of every point, less the start and end 
points, along the route can be computed. With the heading on each point 
established, the change of heading can then be easily computed by taking the 
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difference between the heading of the current point and that of the point prior. 
The computation of the change in heading of the point is illustrated in Figure 57. 



Figure 57. Computation of Change in Heading 

The approach is tested on a short winding route, shown in Figure 58. The 
results of the heading computations are presented in Figures 59 and 60. 







Change in Heading / Deg Heading / Deg 



Figure 59. Heading of Points on Test Route 



Figure 60. Change of Heading of Points on Test Route 
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With the positive results seen from the test, we are able to proceed to 
compute the change in heading for all the points on the actual route. The newly 
computed data for change in heading then allows us to carry out statistical 
analysis to calculate the standard deviation of the distribution. Using the same 
example of the Carmel Valley road, we find the standard deviation of the change 
in heading was calculated to be 5.643 degrees, which classifies the route as a 
type 3 road. This is close to the results obtained by Oriels [10], as discussed in 
Chapter II, and helps to validate the approach adopted in this study. Having 
validated the approach by obtaining similar results with previous works, we then 
proceed to apply the program to routes in other countries where mapping of 
routes may be less extensive than that in the United States. The program is used 
for selected routes in North Korea, Syria, and Afghanistan. The program works 
well, and the results are presented in Figures 61 to Figure 63, respectively. 



Figure 61. Processed North Korean Route—Type 2 
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Figure 62. Processed Syrian Route—Type 2 



Figure 63. Processed Afghanistan Route—Type 3 
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From the results for the routes shown, it can be seen that the approach 
works for the routes that can be selected using Google Maps and imported into 
the software tool. 

The model is initially developed and tested In MATLAB. As discussed 
earlier, a fresh copy is translated into the VBA environment to suit the intended 
user. With the click of a single macro button by the user, the VBA program will 
automatically extract the route data from the GPX file and analyze the route to 
compute the standard deviation in the change of heading, thereby categorizing 
the route into the appropriate road type based on the results obtained. The user 
interface in the VBA environment is presented in Figure 64. 
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Figure 64. Route Analysis in VBA Environment 


Both models are able to produce the same results, verifying that the 
working approach is successfully implemented in an interface suitable for the 
Intended users. 

A copy of the code for the route analysis implemented for MATLAB and 
VBA can be found in Appendix F and Appendix G, respectively. 
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VII. CONCLUSION AND RECOMMENDATIONS 


A. CONCLUSION 

This thesis has achieved the objectives set by successfully generating the 
TSPI files for different synthesized routes that are representative of actual target 
routes of interest. With these TSPI files obtained, the routes can then be fed into 
the weapon IPS model for weapon assessment against a moving target traveling 
along the specific route type. The results of the weapon assessment can then be 
compiled to create a quick look-up reference for assessment of the effectiveness 
of a specific weapon against moving targets moving along route of a particular 
road type. 

A program, written in the VBA environment, is also developed to enable 
quick categorization of an actual road of interest into one of the road profiles 
described, and this program enables the weaponeer to obtain quick estimates of 
the effectiveness of various weapon types available and facilitates a quick 
selection of weapon type for optimal probability of kill. 

The key products of this thesis include: 

• the MATLAB program for generation of synthesized constrained 
routes; 

• the MATLAB program for generating TSPI files for synthesized 
routes.; 

• the Simulink program for generating expected unconstrained 
routes for sea vessels; 

• the MATLAB program for generating TSPI files for sea routes; 

• the MATLAB program for digitization of map routes and 
categorization of map routes into specific road types; and 
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• the VBA program for digitization of map routes and categorization 
of map routes into specific road types. 

B. RECOMMENDATIONS 

This thesis can be further extended into the following areas for future 

works: 

• The velocity computation can be done using the Simulink model 
discussed earlier to better capture the effects of vehicle dynamics 
to produce results that are more representative of the actual vehicle 
performance. 

• The velocity computation can be optimized to determine the optimal 
time for a given vehicle to travel a route in minimal time. 

• The TSPI data obtained can be compared to actual vehicle trials to 
compare the theoretical and empirical results to verify the approach 
and assumptions used. 
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APPENDIX A. MATLAB CODE FOR SYNTHESIZED ROUTE 

GENERATION 


clear all 
clc 

heading_limit=2; % # of sigmas heading is 

limited 

%-New data set for road types 1-6- 

mean_heading= [135 135 135 135 135 135]; 
sigma_heading= [25 29 41 53 66 73]; 
mean_dheading= [0 0 0 0 0 0]; 
sigma_dheading= [0.5 2.5 7 12 17 25]; 
total_km = [40 40 40 40 40 40]; 

road_type=menu( ' Choose a road type 'Type 1', ' Type 

2 ' 

^ / • • • • 

' Type 3 ' , ' Type 4 ' , 'Type 5' , ' Type 6 ' ); 
r_type=[ 'Type 1' ; 'Type 2' ; 'Type 3' ; 'Type 4' ; 'Type 
5' ; 'Type 6' ]; 

r_type_data=cellstr(r_type); 

% get statistics for the chosen road type: 

av_heading = 90 - mean_heading(road_type); % rotate 

coordinates from north to east 

sig_heading=sigma_heading(road_type); 

dhmean = mean_dheading(road_type); 

dhstd_dev = sigma_dheading(road_type); 

h_limit = heading_limit*sig_heading; % heading 

constraint limit 


heading(1) = av_heading; % starting heading 

seg_dist = 10; % 10 meter segments 

npts = total_km(road_type) * 1000 / seg_dist; 


% Randomly sample from the normal distribution 
% Each r value gives you the change in heading for the 
next 10 m of road 
for 1=1:npts 

%r(i) = random('normdhmean,dhstd_dev); 
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r(i)=randn*dhstd_dev+dhmean; 
end 

stdev_data=std(r) % Calculate the actual sigma of 

this data set 


o 

% Here's the road building part, which turns a heading 
change into x and y 

% coordinates starting at easting(l), northing{l). 
easting(l) = 615788; 
northing(l) = 4035674; 

for i=2:npts+1 


heading(i) = heading(i-l) + r(i-l); 


%try to keep road from maki 
if heading(i) > (av_heading 

(av_heading - h_limit) 

heading(i) = heading(i-1) 
end 

easting(i) = easting(i-l) + 
1)*pi/180); 

northing(i) = northing(i-1) 
1)*pi/180); 


ng circles: 

+ h_limit) I heading(i) < 
- r (i-1); 

seg_dist*cos(heading(i- 
+ seg_dist*sin(heading(i- 


end 


% plot the road 
figure(2) 

plot (easting/10 0 0, northing/1 000,' LineWidth ',1.25) 

axis (' equal ') ; grid 

% title('Synthesized Road') 

title(r_type_data(road_type)) 

xlabel( 'Easting (km)' ) 

ylabel( 'Northing (km)' ) 
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APPENDIX B. SIMULINK MODEL FOR VELOCITY COMPUTATION 



77 































































THIS PAGE INTENTIONALLY LEFT BLANK 


78 



APPENDIX C. SIMULINK MODEL FOR SEA ROUTE GENERATION 



Sutisystefn 
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APPENDIX D. MATLAB CODE FOR SYNTHESIZED ROUTE 
PROCESSING INTO TSPI FORMAT 


clear all 
clc 

load( 'Types.mat' ) 
for 1=1:length(types) 

newroute(i,1) = typeS(i,l) - typeS(l,l); % initialise 
1st data point as origin 

newroute(i,2) = typeS(i,2) - typeS(l,2); 
end 

%% Initialization for Computing Speed 
veh_limit =60; % vehicle speed limit in mph 
m2ft = S.28084; % conversion from metres to feet 
mph2ms = 1.609S4 * 1000 / S600; % conversion from mph 
to m/s 

coeff_friction =0.7; % coefficient of friction 

safety_factor = 1; % speed safety factor 

interval = 10; % distance between subsequent points is 

10m 

check =10; % max number of points to check prior and 
after 


decel_limit = 
to m/s2 

0.6* 

9.81; 

% deceleration 

limit, 

from 

g 

accel_limit = 

0.3 * 

9.81; 

% acceleration 

limit, 

from 

g 


to m/s2 

newroute(1,S) = 0; % initialize distance 
newroute(1,4:9) = 0; 

newroute(1,10) = safety_factor * min(sqrt(15 * lel5 * 
m2ft) * (coeff_friction + 0.01 * 0),veh_limit) * 

mph2ms; % initialize vehicle at speed limit, with bank 
angle = 0 

for i=2:length(newroute(:,1))-1 

newroute (i, S) = ( (newroute (i, 1)-newroute (i-1, 1)).'^2 + 
(newroute (i, 2)-newroute (i-1,2) ).''2) .0.5; % there are 

values less than zero as previous computation was based 
on accumulated distance travelled 
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newroute(i,4) = (newroute(i,2) - newroute(i-1,2)) / 

(newroute(i,1) - newroute(i-1,1)); % gradient between 

point i and i-1 

newroute(i,5) = (newroute(i + 1,2) - newroute (i,2)) / 

(newroute(i+1,1) - newroute(i,1)); % gradient between 

point i and i+1 

newroute(i,6) = (newroute(i,4) * newroute(i,5) * 

(newroute(i+1,2) - newroute(i-1,2)) + newroute(i,4) * 

(newroute(i,1) + newroute(i+1,1)) - newroute(i,5) * 
(newroute(i,1) + newroute(i-1,1))) / (2 * 

(newroute(i,4) - newroute(i,5))); % compute x- 

coordinate of centre of circle of curvature 

newroute(i,7) = - 1 / newroute(i,5) * (newroute(i,6) 

- (newroute(i,1) + newroute(i+1,1)) / 2) + 

(newroute(i,2) + newroute(i+1,2)) / 2; % compute y- 
coordinate of centre of circle of curvature 

newroute (i,8) = min(((newroute(i,6) - newroute(i,1)) 

^ 2 + (newroute (i, 7 ) - newroute (i, 2 ) ) 2 ) 0.5,lel5); 

% radius of curvature in metres, lel5 for those with 
infinite values, i.e collinear points 
newroute(i,9) = 0; % road bank angle 
newroute(i,10) = safety_factor * min(sqrt(15 * 
(newroute(i,8) * m2ft) * (coeff_friction + 0.01 * 

newroute(i,9))),veh_limit) * mph2ms; % max speed limit 

at point i in m/s, converted from mph using radius of 
curvature in feet 

newroute(i,16) = sum(newroute(1:i,3)); % Displacement 
end 

%% Generate Speeds at points 

newroute (1, 11) = sqrt (newroute (2,10 )'^2 + 2 * 
decel_limit * interval); 

newroute (1, 12) = sqrt (newroute (2,10) "^2 - 2 * 
accel_limit * interval); 

newroute(1,13) = max(newroute(1,11),newroute(1,12)); 

for i=2:length(newroute(:,!))-! 

% check for deceleration constraint 
acc_checkl = NaN(check,1); 
vel_checkl = NaN (check,1); 

for j=l:min((length(newroute)-1) - i, check) % 
iteration for check on speed limits on points ahead 
if newroute(i,10) > newroute(i+j,10) 
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vel_checkl ( j ) = sqrt (newroute (i +j , 10 )'^2 + 2 * j * 
decel_limit * interval); 
else 
end 
end 

newroute(i,11) = min(min(vel_checkl),newroute(1,10)); 

% check for acceleration constraint 
acc_check2 = NaN(check,1); 
vel_check2 = NaN (check,1); 

for j=l:min(i-l, check) % iteration for check on 
speed limits on points prior 

if newroute(i,10) > newroute(i-j,10) 

vel_check2 ( j ) = sqrt (newroute (i-j , 10 ) "^2 + 2 * j * 
accel_limlt * interval); 
else 
end 
end 

newroute (i,12) = min(min(vel_check2),newroute(1,10)); 
newroute(i,13) = min(newroute(1,11),newroute(i,12)); 

% Actual speed/ highest speed possible 
end 

%% Generate time between points 
for i = l:length(newroute(: , 1))-2 

newroute (i, 14 ) = (newroute (i + 1, 13) "^2 - 
newroute (1, 13)'^2 ) / (2 * interval); % Acceleration 

newroute(i,15) = 2 * interval / (newroute(i,13) + 
newroute(i+1,13)); % time from point i and i+1 
end 

%% Generate Points at 1 sec Intervals 
carryover = 0; 
k = 2; 

time = newroute (:,15); 
newroutetIme(1,1) = newroute(1,1); 
newroutetIme(1,2) = newroute(l,2); 
newroutetlme(1,4) = 0 ; 
for j=l:length(time)-1 

if time(j) + carryover < 1 

carryover = time (j) + carryover; 

else 
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for count=l:floor((carryover + time(j))/I) 
newroutetime(k,8) = k - 1; % total time 
newroutetime(k,1)=newroute(j,1)+(newroute(j+1,1)- 
newroute(j,1))*(count-carryover)/time(j); 

newroutetime(k,2)=newroute(j,2)+(newroute(j+1,2)- 
newroute(j,2))*(count-carryover)/time(j); 

newroutetime(k,9)=newroute(j,13)+(newroute(j+1,13)- 
newroute(j,13))*(count-carryover)/time(j); % speed at 
point k 

k=k+l; 

end 

carryover = carryover + time(j) - count; 
end 
end 

%% Heading and Speed for Time-Based route 
for i=2:length(newroutetime)-1 % heading 
newroutetime(i,4) = newroutetime(i+1,1) - 
newroutetime(i-1,1); % delta x 

newroutetime(i,5) = newroutetime(i+1,2) - 
newroutetime(i-1,2); % delta y 
newroutetime(i,6) = 

atan2(newroutetime(i,5),newroutetime(i,4)); % heading 
with horizontal as 0 rad 

if newroutetime(i,6) > pi/2 && newroutedist(i,6) < pi 
% heading with vertical as 0 rad 

newroutetime(i,7) = 2 * pi + ( pi / 2 - 

newroutetime(i,6)); 
else 

newroutetime(i,7) = pi / 2 - newroutetime(i,6); 
end 

newroutetime(i,10) = newroutetime(i,9) * 

sin(newroutetime(i,7)); 

newroutetime(i,11) = newroutetime(i,9) * 

cos(newroutetime(i,7)); 
end 

%% Conversion From ENU to ECEF to Geodetic 

a = 6378137.0; % WGS-84 Earth semimajor axis (m) 

b = 6356752.3142; % WGS-84 Earth semiminor axis (m) 

f = (a-b)/a; % Ellipsoid Flatness 

e_sq = f * (2 - f); % Square of Eccentricity 
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epsilon = e_sq / (1 - e_sq); 

h = 0; 

routel(l,l) = 0; % initialise synthesized route start 
point lat/lon to 0,0 
routel(1,2) = 0; 

routel(l,3) = sin (routel(1,1)); % sin (lat) 

routel (1,4) = a / sqrt(l-e_sq * routel(1,3) * 

routel(1,3)); % N, distance from the surface to the z- 

axis along the ellipsoid normal 

routel(l,5) = (h + routel(l,4)) * cos (routel (1,1)) * 

cos (routel(1,2)); % converting to ECEF x-coordinates 
routel(1,6) = (h + routel (1,4)) * cos (routel (1,1)) * 

sin(routel(1,2)); % converting to ECEF y-coordinates 
routel (1,7) = (h + (1 - e_sq) * routel (1,4)) * 

sin(routel(1,1)); % converting to ECEF z-coordinates 
for i=l:length(newroutetime) % converting points for Is 
intervals back to Lat/Lon 

route3(i,18) = -sin (routel(1,2)) * newroutetime(i,1) 

- sin(routel (1,1)) * cos (routel (1,2)) * 

newroutetime(i,2) + cos(routel(1,1)) * cos (routel (1,2)) 

* newroutetime(i,3) + routel (1,5); % convert back to 
ECEF x-coordinate 

route3(i,19) = cos(routel(1,2)) * newroutetime(i,1) - 
sin(routel (1,1)) * sin(routel(1,2)) * newroutetime(i,2) 

+ cos (routel (1,1)) * sin (routel (1,2)) * 

newroutetime(i,3) + routel(1,6); % convert back to ECEF 
y-coordinate 

route3(i,20) = cos (routel(1,1)) * newroutetime(i,2) + 
sin (routel (1,1)) * newroutetime(i,3) + routel (1,7); % 

convert back to ECEF z-coordinate 

route3(i,21) = sqrt (route3 (i, 18)'^2 + route3 (i, 19)'^2 ) ; 

% p 

route3(i,22) = atan ( (route3(i,20) * a)/ (route3(i,21) 

* b) ) ; % q 

route3(i,23) = atan2 ( (route3(i,20) + epsilon * b * 
(sin (route3 (i, 22) ) )'^3) , (route3 (i, 21) - e_sq * a * 

(cos (route3 (i, 22) ) ) "^3) ) ; % Lat 

route3(i,24) = atan2 (route3(i,19),route3 (i,18)); % 
Lon 
end 
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NRTLatLon(:,1)=route3(:,23)*180/pi; %converting Is x- 
coordinates back to Lon 

NRTLatLon(:,2)=route3(:,24)*180/pi; %converting Is y- 
coordinates back to Lat 
%% Extracting TSPI 

TSPI = [newroutetime(:,8) NRTLatLon(:,1) NRTLatLon(:,2) 
newroutetime(:,10) newroutetime(:,11)]; 
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APPENDIX E. MATLAB CODE FOR SEA ROUTE PROCESSING 


clear all 
clc 

%% Simulate the Simulink Model 

sim( 'Sea_Vessel_Plant.six' ,90); % runnning Model in 
Simulink 

%% Read the timeseries object from the data file 
load boat.mat 
%% Farcing 

time_sim=boat_results.time; 
route(:,1)=boat_results.Data(:,1); 
route(:,2)=boat_results.Data(:,2); 
rudder = boat_results.Data(:,3); 
speed = boat_results.Data(:,4); 

%% Convert Intervals of 10m 

for 1=1:length(route(:,1))-1 %finding difference in x 

and y for points i and 1+1 

difflx(1)=route (i + 1,1)-route (1,1); 

diffly(1)=route(i+1, 2) -route(1,2); 

end 

k=2; %initiatise k=2 

carryover(1)=0; % initiatise carryover=0 
addx=0; addy=0; 

newroute=[route(1,1) route (1,2)]; 
for j=l:i 

distance (j) = (difflx(j) '^2+dif f ly ( j ) "^2 ) '^0.5; % computing 
distance between points j and j+1 
if carryover(j) > 0 

if (carryover(j)+distance(j))>=10 

for count=1:floor((carryover(j)+distance(j))/lO) 
if count==l 

newroute(k,1)=newroute(k- 

1.1) +difflx(j)/distance(j)*(10-carryover(j))+addx; 

newroute(k,2)=newroute(k- 

1.2) +diffly(j)/distance(j)*(10-carryover(j))+addy; 

k=k+l; 

else 


87 


newroute(k,1)=newroute(k- 

1.1) +difflx(j)/distance(j)*10 ; 

newroute(k,2)=newroute(k- 

1.2) +diffly(j)/distance(j)*10; 

k=k+l; 

end 

end 

carryover(j+1)=distance(j)-(10-carryover(j))- 
(count-1)*10; 

addx=difflx(j)/distance(j)*carryover (j + 1); 
addy=diffly(j)/distance(j)*carryover(j+1); 
else 

carryover(j+1)=carryover(j)+distance(j); 
addx=addx+difflx(j); 
addy=addy+diffly(j); 
end 
else 

if distance(j)>=10 % check if distance to next 
point is > 10m 

for count=1:floor(distance(j)/lO) 
newroute(k,1)=newroute(k- 

1.1) +difflx(j)/distance(j)*10; 

newroute(k,2)=newroute(k- 

1.2) +diffly(j)/distance(j)*10 ; 

k=k+l; 

end 

carryover(j+1)=distance(j)-(count)*10; 
addx=difflx(j)/distance(j)*carryover(j+1); 
addy=diffly(j)/distance(j)*carryover(j+1); 
elseif distance(j)<10 % check if distance to next 
point is < 10m 

carryover(j+1)=distance(j); 
addx=difflx(j); 
addy=diffly(j); 
end 
end 
end 

%% Speed 

for i = l:length(newroute( : , 1) ) 

row_index = find(newroute(i,1)<route(:,1),1)-1; 
if row_index < length(newroute(:,1)) 
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newroute (i,3) = (newroute(i,1) - route(row_index,1)) 

/ (route(row_index+l,1) - route(row_index,1)) * 

(speed(row_index + 1) - speed(row_index)) + 

speed(row_index); % interpolate speed 
end 
end 

%% Generate time between points 
interval = 10; 

for i=l:length(newroute(:,1))-1 

newroute(i,4) = 2 * interval / (newroute(i,3) + 

newroute(i+1,3)); % time from point i and i+1 
end 

%% Generate Points at 1 sec Intervals 
carryover = 0; 
k = 2; 

time = newroute(:,4); 
newroutetime(1,1) = time_sim(l); 
newroutetime(1,2) = newroute(1,1); 
newroutetime(1,3) = newroute(1,2); 
for j=l:length(time)-1 

if time(j) + carryover < 1 

carryover = time (j) + carryover; 

else 

for count=l:floor((carryover + time(j))/I) 
newroutetime(k,1) = k - 1; 

newroutetime(k,2)=newroute(j,1)+(newroute(j+1,1) 
newroute(j,1))*(count-carryover)/time(j); 

newroutetime(k,3)=newroute(j , 2) + (newroute(j + 1,2 ) 
newroute(j,2))*(count-carryover)/time(j); 
k=k+l; 
end 

carryover = carryover + time(j) - count; 
end 
end 

%% Speed 

for i=l:length(newroutetime(:,1)) 

row_index = find(newroutetime(i,2)<route(:,1),1)-1; 
if row_index < length(route (:,1)) 
newroutetime(i,4) = (newroutetime(i,2) - 

route(row_index,1)) / (route(row_index+l,1) - 

route(row_index,1)) * (speed(row_index +1) - 
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interpolate 


speed(row_index)) + speed(row_index); % 

speed 
end 
end 


%% Angle 

newroutetime(1,8) = pi / 2; 

newroutetime(1,9) = newroutetime(1,7) /pi * 180; 
for i=2:length(newroutetime(: , 1))-1 

newroutetime(i,5) = newroutetime(i+1,2) - 
newroutetime(i-1,2); % delta x 

newroutetime(i,6) = newroutetime(i+1,3) - 
newroutetime(i-1,3) ; % delta y 
newroutetime(i,7) = 

atan2(newroutetime(i,6),newroutetime(i,5)); % heading 
with horizontal as 0 rad 

if newroutetime(i,7) > pi/2 && newroutetime(i,7) < pi 

% heading with vertical as 0 rad 

newroutetime(i,8) = 2 * pi + ( pi / 2 - 
newroutetime(i,7)); 
else 

newroutetime(i,8) = pi / 2 - newroutetime(i,7); 
end 

newroutetime(i,9) = newroutetime(i,8) /pi * 180; 
end 

%% Rudder 

for i=l:length(newroutetime(:,1)) 

row_index = find(newroutetime(i,2)<route(:,1),1)-1; 
if row_index < length(route(:,1)) 
newroutetime(i,10) = (newroutetime(i,2) - 
route(row_index,1)) / (route(row_index+l,1) - 

route(row_index,1)) * (rudder(row_index + 1) - 

rudder(row_index)) + rudder(row_index); % interpolate 

speed 
end 
end 

for i=l:length(newroute(:,1)) 

row_index = find(newroute(i,1)<route(:,1),1)-1; 
if row_index < length(route(:,1)) 
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newroute (i,5) = (newroute(i,1) - route(row_index,1)) 

/ (route(row_index+l,1) - route(row_index,1)) * 

(rudder(row_index + 1) - rudder(row_index)) + 

rudder(row_index); % interpolate speed 
end 
end 

TSPI = [newroutetime(:,1) newroutetime(:,2) 
newroutetime(:,3) newroutetime(:,4).* 
sin(newroutetime(:,8)) newroutetime(:,4).* 
cos(newroutetime(:,8))]; 
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APPENDIX F. MATLAB CODE FOR ROUTE ANALYSIS 


clear all 
clc 

%% Generate X Y Z from GPX 

route = gpxread( 'Garmel_Valley.gpx' , 'FeatureType' , 

'track' ); 

%% Defining WGS-84 Geodetic Constants 

a = 6378137.0; % WGS-84 Earth semimajor axis (m) 

b = 6356752.3142; % WGS-84 Earth semiminor axis (m) 

f = (a-b)/a; % Ellipsoid Elatness 

e_sq = f * (2 - f); % Square of Eccentricity 

epsilon = e_sq / (1 - e_sq); 

h = 0; 

routelLat=route.Latitude'; % extracting latitude, phi 
routelLon=route.Longitude'; % extracting longitude, 
lambda 

routelLat2rad=route.Latitude(pi/180); % extracting 
latitude, phi 

routelLon2rad=route.Longitude(pi/180); % extracting 
longitude, lambda 

routel=[routelLat2rad routelLon2rad]; 
for i=l:length(routel) 

routel(i,3) = sin (routel(i,1)); % sin (lat) 
routel(i,4) = a / sqrt(l-e_sq * routel (i,3) * 

routel(i,3)); % N, distance from the surface to the z- 
axis along the ellipsoid normal 

routel (i,5) = (h + routel(i,4)) * cos (routel (i,1)) * 

cos (routel (i,2)); % converting to ECEF x-coordinates 
routel (i,6) = (h + routel(i,4)) * cos (routel (i,1)) * 

sin (routel(i,2)); % converting to ECEF y-coordinates 
routel (i,7) = (h + (1 - e_sq) * routel (i,4)) * 

sin (routel (i,1)); % converting to ECEF z-coordinates 
routel (i,8) = routel(i,5) - routel(1,5); % delta x 
routel(i,9) = routel(i,6) - routel(1,6); % delta y 
routel(i,10) = routel(i,7) - routel(1,7); % delta z 
routel(i,ll) = -sin (routel (1,2)) * routel(i,8) + 

cos (routel (1,2)) * routel(i,9); % Easting 
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routel(i,12) = -cos (routel (1,2)) * sin(routel(1,1)) * 

routel(i,8) - sin(routel(1,1)) * sin(routel (1,2)) * 

routel(i,9) + cos (routel(1,1)) * routel (i,10); % 
Northing 

routel (i,13) = cos (routel(1,1)) * cos (routel (1,2)) * 
routel(i,8) + cos (routel(1,1)) * sin (routel (1,2)) * 

routel(i,9) + sin(routel(1,1)) * routel(i,10); % Up 
end 

%% Convert Intervals of 10m 

for i=l:find(routel(:,11)==routel(end,11))-1 %finding 

difference in x and y for points i and i+1 

dif flx (i)=routel (i + 1,11)-routel (i,11); 

diffly(i)=routel(i+1,12)-routel(i,12); 

end 

k=2; %initiatise k=2 

carryover(1)=0; % initiatise carryover=0 
addx=0; addy=0; 

newroute=[routel(1,11) routel (1,12)]; 
for j=l:i 

distance (j) = (difflx(j) '^2+dif f ly ( j ) '^2 ) ''0.5; % computing 
distance between points j and j+1 
if carryover(j) > 0 

if (carryover(j)+distance(j))>=10 

for count=1:floor((carryover(j)+distance(j))/lO) 
if count==l 

newroute(k,1)=newroute(k- 

1.1) +difflx(j)/distance(j)*(lO-carryover(j))+addx; 

newroute(k,2)=newroute(k- 

1.2) +diffly(j)/distance(j)*(lO-carryover(j))+addy; 

k=k+l; 

else 

newroute(k,1)=newroute(k- 

1.1) +difflx(j)/distance(j)*10 ; 

newroute(k,2)=newroute(k- 

1.2) +diffly(j)/distance(j)*10; 

k=k+l; 

end 

end 

carryover(j+1)=distance(j)-(10-carryover(j))- 
(count-1)*10; 

addx=difflx(j)/distance(j)*carryover(j+1); 
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addy=diffly(j)/distance(j)*carryover(j+1); 
else 

carryover(j+1)=carryover(j)+distance(j); 
addx=addx+difflx(j); 
addy=addy+diffly(j); 
end 
else 

if distance(j)>=10 % check if distance to next 
point is > 10m 

for count=1:floor(distance(j)/lO) 
newroute(k,1)=newroute(k- 

1.1) +difflx(j)/distance(j)*10; 

newroute(k,2)=newroute(k- 

1.2) +diffly(j)/distance(j)*10 ; 

k=k+l; 

end 

carryover(j + 1)=distance(j)-(count)*10 ; 
addx=difflx(j)/distance(j)*carryover(j+1); 
addy=diffly(j)/distance(j)*carryover(j+1); 
elseif distance(j)<10 % check if distance to next 
point is < 10m 

carryover(j+1)=distance(j); 
addx=difflx(j); 
addy=diffly(j); 
end 
end 
end 


%% Computing UP Coordinates from Northing 
rindex=0; 

newroutedist(:,1) = newroute(:,1); 
newroutedist(:,2) = newroute(:,2); 

for i=l:length(newroutedist) % Computing UP Coordinates 

from Northing for 10m interval 

rindex=find(newroutedist(i,2)>=routel(:,12),1); 
newroutedist(i,3)=(newroutedist(i,2)- 
routel(rindex,12))/(routel(rindex+1,12)- 
routel(rindex,12)) * (routel(rindex+1,13)- 

routel(rindex,13)) + routel(rindex,13); 
end 
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%% Heading Computation for Distance-Based 

for i=2:length(newroutedist)-1 

newroutedist(i,4) = newroutedist(i+1,1) - 
newroutedist(i-1,1); % delta x 

newroutedist(i,5) = newroutedist(i+1,2) - 
newroutedist(i-1,2); % delta y 
newroutedist(i,6) = 

atan2(newroutedist(i,5),newroutedist(i,4)); % heading 
with horizontal as 0 rad 

if newroutedist(i,6) > pi/2 && newroutedist(i,6) < pi 
% heading with vertical as 0 rad 

newroutedist(i,7) = 2 * pi + ( pi / 2 - 

newroutedist(i,6)); 
else 

newroutedist(i,7) = pi / 2 - newroutedist(i,6); 
end 
end 

for i=2:length(newroutedist)-2 

if newroutedist(i+1,7) - newroutedist(i,7) >= pi && 

newroutedist(i+1,7) - newroutedist(i,7) <= 3 * pi % 

change in heading 

newroutedist(i,8) = newroutedist(i+1,7) - 

newroutedist(i,7) - 2 * pi; 

elseif newroutedist(i+1,7) - newroutedist(i,7) <= -pi 

&& newroutedist(i+1,7) - newroutedist(i,7) >= -3 * pi 

newroutedist(i,8) = newroutedist(i+1,7) - 

newroutedist(i,7) + 2 * pi; 

else 

newroutedist(i,8) = newroutedist(i+1,7) - 

newroutedist(i,7); 

end 

end 

if std(newroutedist(2:length(newroutedist)-2,8) )/pi*180 
>= 0 && std(newroutedist(2:length(newroutedist)- 
2,8))/pi*180 <= 1 % classify road type 
roadtype = 1; 

elseif std(newroutedist(2:length(newroutedist)- 
2,8))/pi*180 > 1 && 

std(newroutedist(2:length(newroutedist)-2,8))/pi*180 <= 
4.5 
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roadtype = 2; 

elseif std(newroutedist(2:length(newroutedist)- 

2,8))/pi*180 > 4.5 && 

std(newroutedist(2:length(newroutedist)-2,8))/p1*180 < 

9.5 

roadtype = 3; 

elseif std(newroutedist(2:length(newroutedist)- 

2,8))/pi*180 > 9.5 && 

std(newroutedist(2:length(newroutedist)-2,8))/p1*180 < 

14.5 

roadtype = 4; 

elseif std(newroutedist(2:length(newroutedist)- 

2,8))/pi*180 > 14.5 && 

std(newroutedist(2:length(newroutedist)-2,8))/p1*180 < 

19.5 

roadtype = 5; 

elseif std(newroutedist(2:length(newroutedist)- 

2,8))/pi*180 > 19.5 && 

std(newroutedist(2:length(newroutedist)-2,8))/p1*180 > 

30.5 

roadtype = 6; 
else 

roadtype = 7; 
end 

fprintf( '\nThe sigma change in heading is %g degrees, 
therefore it is road type %d.\n', 

std(newroutedist(2:length(newroutedist)-2,8))/p1*180, 
roadtype) 
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APPENDIX G. VBA CODE FOR ROUTE ANALYSIS 


Sub Main() 

With ThisWorkbook 
■SheetsC'Workspace"). Visible = True 
Call ImportGPX 
Call CopyColumnToWorkbook 
Call FillDownXY 
Call convertm 
Call HeadingXY 
Call Chart 
Call RemoveGPX 

'.SheetsC'Workspace").Visible = xlSheetVeryHidden 
End With 
End Sub 

Sub ImportGPXO ' Import GPX data into worksheet 
With ThisWorkbook 

.Sheets.Add(After:=.Sheets(.Sheets.Count)).Name = "GPX data" 

End With 

Application.DisplayAlerts = False 

ActiveWorkbook.Xml Import Url:="C:\Weapon\GPXfile\Route.gpx", 

lmportMap:=Nothing, _ 

Overwrite:=True, Destination:=Worksheets("GPX data").Range("$A$1") 
Application.DisplayAlerts = True 

End Sub 

Sub CopyColumnToWorkbookO 
Dim r1 As Range 
With ThisWorkbook 
.SheetsC'GPX data").Select 
Set r1 = Range("L:M") 

.SheetsC'Workspace") .Select 

.SheetsC'Workspace").Range("A4:Z1000000").ClearContents 
rl.Copy Range("A1") 

End With 
End Sub 


Sub FillDownXYO 

Dim strFormulasI (1 To 2), strFormulas2(1 To 3) As Variant 
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Dim LastRow As Long 

Dim a, b, f, e_sq, epsilon, h As Long 

a = 6378137 'WGS-84 Earth semimajor axis (m) 

b = 6356752.3142 'WGS-84 Earth semiminor axis (m) 

f = (a - b) / a 'Ellipsoid Flatness 
e_sq = f * (2 - f) 'Square of Eccentricity 
epsilon = e_sq / (1 - e_sq) 


h = 0 


LastRow = Range("A" & Rows.Count).End(xlLlp).Row 

With ThisWorkbook.Sheets("Workspace") 
strFormulasl(l) = "=A4/180*pi()" 
strFormulasI (2) = "=B4/180*pi()" 

.Range("C4:D4").Formula = strFormulasI 
.Range("C4:D" & LastRow).FillDown 

For j = 1 To LastRow - 3 

.RangeC'E" & j -i- 3).Value = a / ((1 - e_sq * Sin(.Range("C" & ] -i- 3).Value) * 
Sin(.Range("C" & j -i- 3).Value)) * 0.5) 'N, distance from the surface to the z-axis 
along the ellipsoid normal 

.RangeC'F" & j -i- 3).Value = (h -i- .Range("E" & j -i- 3).Value) * Cos(.Range("C" 
& j -I- 3).Value) * Cos(.Range("D" & ] -i- 3).Value) 'converting to ECEF x- 
coordi nates 

.RangeC'G" & ] -i- 3).Value = (h -i- .Range("E" & j -i- 3).Value) * Cos(.Range("C" 
& j -I- 3).Value) * Sin(.Range("D" & j -i- 3).Value) 'converting to ECEF y- 
coordi nates 

.RangeC'FI" & ] -i- 3).Value = (h -i- (1 - e_sq) * .Range("E" & j -i- 3).Value) * 
Sin(.Range("C" & ] -i- 3).Value) 'converting to ECEF z-coordinates 

.RangeC'l" & ] -i- 3).Value = .Range("F" & j -i- 3).Value - .Range("F" & 1 -i- 
3).Value 'delta x 

.RangeC'J" & j + 3).Value = .Range("G" & j -r 3).Value - .Range("G" & 1 -r 
3).Value 'delta y 

.RangeC'K" & ] -i- 3).Value = .Range("H" & ] -i- 3).Value - .Range("FI" & 1 -i- 
3).Value 'delta z 

.RangeC'L" & ] + 3).Value = -Sin(.Range("D" & 1 + 3).Value) * .Range("l" & ] + 
3).Value -i- Cos(.Range("D" & 1 -i- 3).Value) * .Range("J" & j -i- 3).Value 'Easting 

.RangeC'M" & ] -i- 3).Value = -Cos(.Range("D" & 1 -i- 3).Value) * 
Sin(.Range("C" & 1 -i- 3).Value) * .Range("l" & ] -i- 3).Value - Sin(.Range("C" & 1 -i- 
3).Value) * Sin(.Range("D" & 1 -i- 3).Value) * .Range("J" & j -i- 3).Value -i- 
Cos(.Range("C" & 1 -i- 3).Value) * .Range("K" & j -i- 3).Value 'Northing 
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■ RangeC'N" & j + 3).Value = Cos(.Range("C" & 1 + 3).Value) * 
Cos(.Range("D" & 1 + 3).Value) * .Range("l" & j + 3).Value + Cos(.Range("C" & 1 
+ 3).Value) * Sin(.Range("D" & 1 + 3).Value) * .Range("J" & j + 3).Value + 
Sin(.Range("C" & 1 + 3).Value) * .Range("K" & j + 3).Value 'Up 
Next j 

strFormulas2(1) = "=L5-L4" 
strFormulas2(2) = "=M5-M4'' 
strFormulas2(3) = "=SQRT(04'^2+P4^2)" 

.Range("04:Q4").Formula = strFormulas2 
.Range("04:Q" & LastRow - 1).FillDown 
End With 

End Sub 


Sub convertmO 

Dim I, j, k, Count As Integer 

Dim carryover, addx, addy, distance As Long 

addx = 0 
addy = 0 
k = 2 

I = RangeC'P" & Rows.Count).End(xlUp).Row - 3 

With ThisWorkbook.Sheets("Workspace") 

.Range("R4'').Value = 0 
.Range("S4").Value = .Range("L" & 4).Value 
.Range("T4").Value = .Range("M" & 4).Value 
Forj = 1 Toi 

If .RangeC'R" & j + 3).Value > 0 Then 
If .RangeC'R" & j + 3).Value + .Range("Q" & j + 3).Value >= 10 Then 
For Count = 1 To Application.WorksheetFunction.Floor((.Range("R" & j + 
3).Value + .Range("Q" & j + 3).Value) /10, 1) 

If Count = 1 Then 

.RangeC'S" & k + 3).Value = .Range("S" & k + 3 - 1).Value + .Range("0" 
& j + 3).Value / .Range("Q" & j + 3).Value * (10 - .Range("R" & j + 3).Value) + 
addx 

.RangeC'T" & k + 3).Value = .Range("T" & k + 3 - 1).Value + .Range("P" 
& j + 3).Value / .Range("Q" & j + 3).Value * (10 - .Range("R" & j + 3).Value) + 
addy 

k = k+ 1 
Else 
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■ RangeC'S" & k + 3).Value = .Range("S" & k + 3 - 1).Value + .Range("0" 
& j + 3).Value / .Range("Q" & j + 3).Value * 10 

.RangeC'T" & k + 3).Value = .Range("T" & k + 3 - 1).Value + .Range("P" 
& j + 3).Value / .Range("Q" & j + 3).Value * 10 
k = k+ 1 

End If 
Next Count 

■ RangeC'R" & j + 3 + 1).Value = .Range("Q" & j + 3).Value - (10 - 
.RangeC'R" & j + 3).Value) - (Count - 2) * 10 

addx = .RangeC'O" & j + 3).Value / .Range("Q" & j + 3).Value * .Range("R" 
& j + 3 + 1).Value 

addy = .Range("P" & j + 3).Value / .Range("Q" & j + 3).Value * .Range("R" 
& j + 3 + 1).Value 
Else 

■ RangeC'R" & j + 3 + 1).Value = .Range("Q" & j + 3).Value + .Range("R" & j 
+ 3).Value 

addx = .RangeC'O" & j + 3).Value + addx 
addy = .Range("P" & j + 3).Value + addy 
End If 
Else 

If.RangeC'O" & j + 3).Value >= 10 Then 
For Count = 1 To Application.WorksheetFunction.Floor((.Range("Q" & j + 
3).Value)/10, 1) 

■ RangeC'S" & k + 3).Value = .Range("S" & k + 3 - 1).Value + .RangeC'O" 
& j + 3).Value / .RangeC'O" & j + 3).Value * 10 

■ RangeC'T" & k + 3).Value = .Range("T" & k + 3 - 1).Value + .Range("P" & 
j + 3).Value / .RangeC'O" & j + 3).Value * 10 

k = k+ 1 
Next Count 

■ RangeC'R" & j + 3 + 1).Value = .RangeC'O" & j + 3).Value - (Count - 1) * 10 
addx = .RangeC'O" & j + 3).Value / .RangeC'O" & j + 3).Value * .RangeC'R" 

& j + 3 + 1).Value 

addy = .Range("P" & j + 3).Value / .RangeC'O" & j + 3).Value * .RangeC'R" 
& j + 3 + 1).Value 
Else 

If .RangeC'O" & j + 3).Value < 10 Then 

■ RangeC'R" & j + 3 + 1).Value = .RangeC'O" & j + 3).Value 

addx = .RangeC'O" & j + 3).Value 

addy = .Range("P" & j + 3).Value 
End If 
End If 
End If 
Next j 
End With 
End Sub 
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Sub ChartQ 

Dim myChtObj As ChartObject 
Dim rngChtData As Range 
Dim rngChtXVal As Range 
Dim rngRefData As Range 
Dim rngRefXVal As Range 
Dim RngToCover As Range 
Dim LastRow As Long 

With ThisWorkbook 
SheetsC'Workspace"). Activate 
LastRow = RangeC'T" & Rows.Count).End(xlLlp).Row 

Set rngChtData = ActiveSheet.Range("T4:T" & LastRow) 
Set rngChtXVal = ActiveSheet.Range("S4:S" & LastRow) 
Set rngRefData = ActiveSheet.Range("AA4:AA5") 

Set rngRefXVal = ActiveSheet.Range("Z4:Z5") 
Sheets("Main").Activate 

Set myChtObj = ActiveSheet.ChartObjects.Add _ 
(Left:=250, Width:=375, Top:=75, Height:=225) 

With myChtObj.Chart 

.ChartType = xIXYScatterLines 

' remove extra series 
Do Until .SeriesCollection.Count = 0 
.SeriesCollection.Delete 
Loop 

' add series from selected range, column by column 
With .SeriesCollection.NewSeries 
.Values = rngChtData 
.XValues = rngChtXVal 
.Name = "Route Plot" 

End With 

With .SeriesCollection.NewSeries 
.Values = rngRefData 
.XValues = rngRefXVal 
.Name = "1km Reference" 

End With 

Set RngToCover = ActiveSheet.Range("G8:023") 
myChtObj.Height = RngToCover.Height' resize 
myChtObj.Width = RngToCover.Width ' resize 
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myChtObj.Top = RngToCover.Top ' reposition 
myChtObj.Left = RngToCover.Left ' reposition 
End With 
End With 
End Sub 

Sub HeadingXYQ 
Dim I As Integer 

Dim strFormulasI (1 To 3), strFormulas2(1 To 1), strFormulas3(1 To 1), 
strFormulas4(1 To 1), strFormulas5(1 To 1), strFormulas6(1 To 1) As Variant 
Dim LastRow As Long 

With ThisWorkbook 
.SheetsC'Workspace").Select 
End With 

LastRow = RangeC'T" & Rows.Count).End(xlLlp).Row 

With ThisWorkbook.Sheets("Workspace") 

StrFormulasI (1) = "=S6-S4" 

StrFormulasI (2) = "=T6-T4" 

StrFormulasI (3) = "=atan2(LI5,V5)/pi()*180" 

.Range("LI5:W5").Formula = strFormulasI 
.Range("LI5:W" & LastRow - 1).FillDown 
For I = 1 To LastRow - 5 

If .Range('W" & I + 4).Value > 90 And .Range("W" & I + 4).Value < 180 Then 
.RangeC'X" & I + 4).Value = 360 + (90 - .Range("W" & I + 4).Value) 

Else 

.RangeC'X" & I + 4).Value = 90 - .Range("W" & I + 4).Value 
End If 
Next I 

strFormulas2(1) = "=IF(AND(X6-X5<=540,X6-X5>=180),X6-X5- 

360JF(AND(X6-X5>=-540,X6-X5<=-180),X6-X5+360, X6-X5))" 

.Range("Y5").Formula = strFormulas2 
.Range("Y5:Y" & LastRow - 2).FillDown 
.Range("V" & LastRow + 2).Value = "Mean Fleading" 

.RangeC'W" & LastRow + 2).Value = "Sigma Heading" 

.RangeC'X" & LastRow + 2).Value = "Mean Change Heading" 

.RangeC'Y" & LastRow + 2).Value = "Sigma Change Heading" 

.RangeC'Y" & LastRow + 4).Value = "Road Type" 
strFormulas3(1) = "=AVERAGE(X4:X" & LastRow -1 & ")" 
strFormulas4(1) = "=STDEVA(X4:X" & LastRow - 1 & ")" 
strFormulas5(1) = "=AVERAGE(Y5:Y" & LastRow -1 & ")" 
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strFormulas6(1) = "=STDEVA(Y5:Y" & LastRow -1 & ")" 

■ RangeC'V" & LastRow + 3).Formula = strFormulasS 
.Range("W" & LastRow + 3).Formula = strFormulas4 
.Range("X" & LastRow + 3).Formula = strFormulasS 
.Range("Y" & LastRow + 3).Formula = strFormulasS 

If .RangeC'Y" & LastRow + 3).Value >= 0 And .Range("Y" & LastRow + 
3).Value <= 1 Then 

■ RangeC'Y" & LastRow + 5).Value = 1 

Elself .RangeC'Y" & LastRow + 3).Value > 1 And .RangeC'Y" & LastRow + 
3).Value <= 4.5 Then 

■ RangeC'Y" & LastRow + 5).Value = 2 

Elself .RangeC'Y" & LastRow + 3).Value > 4.5 And .RangeC'Y" & LastRow + 
3).Value <= 9.5 Then 

■ RangeC'Y" & LastRow + 5).Value = 3 

Elself .RangeC'Y" & LastRow + 3).Value > 9.5 And .RangeC'Y" & LastRow + 
3).Value <= 14.5 Then 

■ RangeC'Y" & LastRow + 5).Value = 4 

Elself .RangeC'Y" & LastRow + 3).Value > 14.5 And .RangeC'Y" & LastRow + 
3).Value <= 19.5 Then 

■ RangeC'Y" & LastRow + 5).Value = 5 

Elself .RangeC'Y" & LastRow + 3).Value > 19.5 And .RangeC'Y" & LastRow + 
3).Value >= 30.5 Then 

■ RangeC'Y" & LastRow + 5).Value = 6 
Else 

■ RangeC'Y" & LastRow + 5).Value = "ERROR!!" 

End If 

If .RangeC'V" & LastRow + 3).Value >= 0 And .Range("V" & LastRow + 
3).Value <= 90 Then 

.Range("Z4").Value = .Range("S4").Value 
.Range("AA4").Value = .Range("T4").Value - 2000 
.Range("Z5").Value = .Range("S4").Value + 1000 
.Range("AA5").Value = .Range("T4").Value - 2000 

Elself .RangeC'V" & LastRow + 3).Value > 90 And .RangeC'V" & LastRow + 
3).Value <= 180 Then 

.Range("Z4").Value = .Range("S4").Value 
.Range("AA4").Value = .Range("T4").Value + 2000 
.Range("Z5").Value = .Range("S4").Value + 1000 
.Range("AA5").Value = .Range("T4").Value + 2000 

Elself .RangeC'V" & LastRow + 3).Value < 0 And .RangeC'V" & LastRow + 
3).Value >= -90 Then 

.Range("Z4").Value = .Range("S4").Value - 1000 
.Range("AA4").Value = .Range("T4").Value - 2000 
.Range("Z5").Value = .Range("S4").Value 
.Range("AA5").Value = .Range("T4").Value - 2000 


105 



Else 

.Range("Z4").Value = .Range("S4").Value - 1000 
.Range("AA4").Value = .Range("T4").Value + 2000 
.Range("Z5").Value = .Range("S4").Value 
.Range("AA5").Value = .Range("T4").Value + 2000 
End If 
End With 

Worksheets("Workspace").Range("Y" & LastRow + 
Worksheets("Main").Range("F5") 
Worksheets("Main").Range("F5").lnterior.Colorlndex = 6 
Worksheets("Main").Range("F5").Font.Name = "Calibri" 
Worksheets("Main").Range("F5").Font.Size = 20 
Worksheets("Main").Range("F5").Font.Bold = True 
Worksheets("Main").Range("F5").HorizontalAlignment = xlCenter 
End Sub 


Sub RemoveGPXO 
With ThisWorkbook 

Application.DisplayAlerts = False 
.SheetsC'GPX data").Delete 
Application.DisplayAlerts = True 
'.Sheets("Workspace").Visible = xlSheetVeryHidden 
.SheetsC'Main").Select 

End With 
End Sub 


5).Copy 
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